Does the recent revival of Western disturbances govern the Karakoram Anomaly?

: The global retreat in glaciers is considered to be one of the critical indicators of climate change. However, the glaciers of the Karakoram (KR) region of the Karakoram – Himalayas (KH) stand out because of their divergent response, displaying a surge as opposed to glaciers in other regions. This phenomenon is known as the “ Karakoram anomaly. ” Although many factors control the establishment and sustenance of the anomaly, the present study establishes winter precipitation associated with western disturbances (WDs) over the KH as one of the key drivers behind its emergence. To examine the role of WDs, a tracking algorithm is applied to 39 seasons (November – March) for three separate (ERA5, MERRA-2, and NCEP-CFSR/CFSv2) reanalysis datasets. The associated reanalysis ensemble statistics of WD properties produced in terms of their intensity, precipitation/snowfall volumes, and wind speed suggest a revival in recent years over the core-anom-aly regions. However, the frequency has remained steady. The Karakoram has witnessed a rise of ∼ 10% in precipitation intensity associated with WDs in recent decades. The high percentage of snowfall received by the Karakoram ( ∼ 65%) from WDs relative to the total seasonal snowfall suggests a crucial role in modulating the regional mass-balance anomaly. Simulta-neously, the amount of snowfall from non-WD sources in the Karakoram has had a statistically signi ﬁ cant decline of ∼ 17% in recent decades, coinciding with the anomaly period. The enhanced intensity of WDs is found to be associated with changes in increased baroclinic instability and a shift of the subtropical westerly jet mean latitudinal position.


Introduction
Western disturbances (WDs) are upper-tropospheric synoptic-scale cyclonic systems embedded in the subtropical westerly jet stream, often associated with extreme rainfall events in northern India during boreal winter (Hunt et al. 2018b).They are further enhanced over the Karakoram-Himalayas (KH) due to orographic uplift (Hunt et al. 2018b).WD-associated snowfall is the dominant precipitation over the KH during winter, playing a critical role in establishing and sustaining the regional snowpack and, at the same time, replenishing regional water resources (Dimri et al. 2015).Glaciers, one of the critical indicators of global warming, suffer significant impacts due to climate change.The KH region contains the largest glacier mass outside the Arctic and Antarctica.It is often referred to as the "water tower of Asia" (Bolch et al. 2012) and is the source of several perennial rivers that flow through the major basins (Frey et al. 2014) surrounding the "High Asian" region.Even the slightest variation in the annual glacier-melt runoff can have profound impacts on water resources.Understanding the KH glaciers' behavior and their contribution to the sustainable water supply is a significant challenge for the scientific community and policymakers.The KH also significantly impacts the climate and hydrology of one of the world's most densely populated regions.Studying its glaciers is essential for millions of dwellers downstream, especially during the dry season (Immerzeel et al. 2010;Kaser et al. 2010).As a parameter, the temperature usually remains at the forefront when studying changes associated with the glacier mass balance.However, changes in precipitation also play a crucial role in affecting them.
The spatiotemporal distribution of precipitation over the KH is hugely variable and depends on synoptic (and mesoscale) factors and the complex topography (Norris et al. 2015(Norris et al. , 2017(Norris et al. , 2018)).However, their highly fluctuating and erratic precipitation behavior is dominated by two large-scale circulation systems: Indian summer monsoon rainfall (ISMR) and WDs (Bookhagen and Burbank 2010).The ISMR is responsible for most of the precipitation in the central and eastern Himalayas with a relatively subdued impact in the western and Karakoram ranges, which are the primary recipient of most of the precipitation provided by the WDs.Obtaining in situ measurements or setting up weather stations to gauge the impact of WDs over the KH is extremely difficult because of its rugged terrain and topography, limiting the monitoring of its hydrometeorology and microclimate.The hydrological impacts of the glaciers have been widely reported, and they play defining roles in controlling the amount of meltwater contributing to runoff (Immerzeel et al. 2010;M ölg et al. 2014).
Although the global glacial melt is persisting with its positive trend, the Karakoram and Kunlun regions have shown exceptional stability (Forsythe et al. 2017;Kumar et al. 2015;Brun et al. 2017;Hewitt 2005;Gardelle et al. 2012Gardelle et al. , 2013)).This has been identified as the "Karakoram anomaly," for which many hypotheses have been proposed since the early 1990s (Zhou et al. 2017;Minora et al. 2013;Forsythe et al. 2017;Li et al. 2018;Kapnick et al. 2014;de Kok et al. 2018;Farinotti et al. 2020).More profound insights into the mass-balance changes and runoff estimates have been provided by studies such as those on the Chhota Shigri glacier using either climate models (Engelhardt et al. 2017a,b) or various other robust techniques such as surface energy balance measurements or reconstruction of mass balance (Azam et al. 2014a,b).A recent study showed that snowfall variability dictates the mass-balance variability over the KH region (Kumar et al. 2019).Despite ongoing research, questions remain, such as the role of WDs, whose impact is synoptic in nature, in modulating the regional hydrology and glacial anomaly (Norris et al. 2015;Cannon et al. 2016;Norris et al. 2017Norris et al. , 2018)).The requirement is to devise a holistic view of the region, broadly defining its characteristics and societal impacts.A previous study addressed WD dynamics using a wave-tracking algorithm based on 500-hPa geopotential height anomaly and computed around 5 or 6 events per year through northern India (Cannon et al. 2016).However, this number is significantly less than the ones reported elsewhere in the literature.Some studies have suggested the frequency of WDs to be about 4-7 month 21 in boreal winter (Dimri 2006;Midhuna et al. 2020).Another study identified 3090 tracks in 37 years passing through northern India (Hunt et al. 2018b).The differences are likely due to different identification methodologies and different sampling regions.Although WDs can be present at any given time of the year, their occurrence peaks during the winter along with their associated precipitation, mainly consisting of snowfall, which is the most crucial parameter for glacier mass budget estimation during the accumulation period (Kumar et al. 2019).
The aims of this study are 1) to examine the role of WDs in controlling the overall regional glacier mass-balance variability, 2) to understand the influence of WDs behind the emergence and continual presence of the Karakoram anomaly in recent decades, 3) to quantify the contribution of WDs to regional precipitation regimes during the accumulation period, 4) to identify the key behavioral changes in WD frequency and intensity (in terms of precipitation and wind speed), along with their subsequent impact in modulating the eventual mass-balance estimates, and 5) to explore the possible dynamical changes controlling the formation and propagation of WDs in recent decades.The study starts with the description of the study region (section 2a), followed by an outline of the data sources (section 2b) and the methodology employed (section 2c).The results section primarily quantifies the impact of WDs in terms of the changes in their frequency (section 3a), precipitation intensity (section 3b), the volume of WD-associated precipitation and snowfall when compared to total seasonal volumes (section 3c), and the associated wind speed intensity (section 3d).A discussion of recent dynamical changes behind the genesis and distribution of WDs is included in section 3e.The subsequent impact of all the changes on the mass balance of KH glaciers is discussed in section 3f.We summarize and conclude the discussion in section 4.

a. Study region
The KH is divided into four subregions, namely the Karakoram (KR), western Himalayas (WH), central Himalayas (CH), and eastern Himalayas (EH).The boundaries for these regions have been defined in various previous studies (Bolch et al. 2012;Kumar et al. 2015).Figure 1a shows the four regions, along with their topographical features and locations of observation stations.The scarcity of India Meteorological Department (IMD) measuring stations over the highest elevated regions is evident, with most observational stations concentrated at the southern edges of the western and central Himalayas.Similarly, Climate Research Unit (CRU) land stations are almost entirely absent from the Karakoram and central Himalayas.Figure 1b depicts the simulated mass balance of the Karakoram-Himalayan region using a high-resolution dynamically coupled glacier climate model REMO glacier (the details of the model are provided in later sections).The red oval in the figure illustrates how the Karakoram is characterized by anonymously positive mass-balance estimates, diverging from most of the region's glaciers, also shown in Fig. 3 of Bolch et al. (2012) and Fig. 1 of Azam et al. (2018) using field measurements.

b. Datasets
The present study uses three separate reanalysis datasets to identify WDs storm tracks and their associated statistics: ERA5, MERRA-2, and NCEP-CFSR/CFSv2.TRMM precipitation data are used for a comparative assessment of track associated precipitation against the three reanalysis datasets.

1) ERA5
The ERA5 dataset (Hersbach et al. 2020) is the latest reanalysis from the European Centre for Medium-Range Weather Forecasts (ECMWF).This new reanalysis replaces the popular ERA-Interim (Dee et al. 2011) and has a significantly enhanced horizontal resolution (TL639, spectral model, 31 km) and vertical resolution (137 levels to 0.01 hPa, hybrid sigma-pressure) and an hourly temporal output.Data from 1980 to the present are used here.
2) MERRA-2 The Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA-2) (Gelaro et al. 2017), is a global atmospheric reanalysis produced by the NASA Global Modeling and Assimilation Office (GMAO).It spans the satellite observing era from 1980 to the present.It replaces its predecessor, MERRA (Rienecker et al. 2011).
MERRA-2 has a horizontal resolution of 0.6258 3 0.58 and 72 hybrid-eta levels up to 0.01 hPa.Past observations are assimilated using a 3D-Var system with a 6-h update cycle, which only adjusts model projections forward in time.

3) NCEP-CFSR/CFSv2
The third reanalysis used in the study is the NCEP Climate Forecast System Reanalysis (CFSR) (Saha et al. 2010), which provides data from 1979 to 2010.NCEP-CFSR is a coupled reanalysis with a horizontal resolution of ∼38 km (T382, spectral model) and 64 vertical levels (sigma-pressure hybrid) up to 0.3 hPa in the atmosphere.Historical observations were assimilated into the CFS model using a 3D-Variational system, similar to MERRA-2.The upgraded successor of CFSR, namely the NCEP Climate Forecast System Reanalysis version 2 (CFSv2) (Saha et al. 2014), extends it from 2011 to the present.CFSR can be accessed from  Using three separate reanalyses helps quantify and explain the associated uncertainties involved in identifying and tracking the WDs and their associated statistics.To achieve the objectives, the present work utilizes 3-hourly global relative vorticity fields at 400 and 300 hPa (Hunt et al. 2018a(Hunt et al. , 2019a,b) ,b) for the common period  to identify and track the WDs as well as to examine their properties.The two levels are vertically averaged to get the best signals at approximately 350 hPa (Hunt et al. 2019a).Other parameters used in the study (either individually or as the ensemble mean) are the total precipitation, snowfall, dewpoint temperature, u and y winds, and atmospheric temperature at 400 and 300 hPa.

c. Methodology
The present study aims to quantify the impact of WDs on the Karakoram-Himalayan glaciers in terms of their frequency and intensity in recent years.The study period  has been divided into two subperiods: P1 (1980-2000) and P2 (2001-19).The choice of periods is inspired by the observed glacier variability in the region.Various studies have discussed the remarkable change in the trends of mass balance of Karakoram-Himalayan glaciers at the turn of the century (Gardelle et al. 2012;Bolch et al. 2012;Dehecq et al. 2019;Berthier and Brun 2019).By comparing and analyzing various parameters within the two subperiods, the study attempts to unravel the reasons behind this behavioral change of KH glaciers since the start of the twenty-first century, including the much-debated Karakoram anomaly.

1) TRACKING ALGORITHM
A tracking algorithm was applied to the 3-hourly uppertropospheric vorticity field, averaged between 400 and 300 hPa, for 39 winter seasons in the three reanalyses to understand the impact and contribution of WDs to the regional snow and ice replenishment.The tracking is performed for the November to March season, considered the most active period for WDs (Dimri et al. 2015;Krishnan et al. 2018).The tracking algorithm ("TRACK,") has been extensively used to study various weather systems across the globe (Hoskins and Hodges 2002;Hodges et al. 2017;Curio et al. 2019;Pinheiro et al. 2019;Rastogi et al. 2018).The pressure levels used in the present study were carefully chosen after analyzing the spectrum of levels (between 500 and 200 hPa) used in previous studies (Hunt et al. 2018b;Krishnan et al. 2018;Madhura et al. 2014) to explore the dynamics of these synoptic-scale phenomena, keeping in mind that the typical WDs have their peak vorticity at 350 hPa (Hunt et al. 2019b).Vorticity is preferred over other identification variables such as geopotential as it focuses on smaller spatial scales.
Vorticity is directly available from ERA5 as it is a prognostic variable in that model; however, for the other reanalysis where only winds are available, the relative vorticity is easily computed as the vertical component of the curl of the winds in spherical coordinates: where u is latitude, k is longitude, u is the zonal wind, and y is the meridional wind.This is performed using B-splines from which the first-order derivatives can be obtained analytically rather than using finite differences.Since vorticity is a noisy field at the native resolutions of the reanalyses, spectral filtering was applied to focus on synoptic-scale features.The vertically averaged vorticity is spatially filtered using the spectral filtering to T63 resolution (triangular truncation at total wavenumber 63, isotropic).The large-scale background is also removed for total wavenumbers less than 5. Additionally, the spectral coefficients are tapered to smooth the data further and suppress Gibbs oscillations (Sardeshmukh and Hoskins 1984).The vorticity maxima (cyclones) are initially determined on the T63 grid that exceeds a threshold value of 1.0 3 10 25 s 21 and then used as starting points to obtain the off-grid locations using B-spline interpolation and maximization methods (Hodges 1995), resulting in smoother tracks.Tracks are first initialized using the identified maxima based on the nearest neighbor method and then refined by minimizing a cost function for track smoothness subject to adaptive constraints for displacement and track smoothness (Hodges 1995(Hodges , 1999) ) with the constraints chosen appropriately for the 3-hourly data.This results in a coherent set of tracks for full system life cycles, including their nascent, mature, and decay stages.After the tracking, the corresponding full-resolution total precipitation, snowfall, and vertically averaged wind speed at 400-300-hPa levels are associated with each point of the tracks for different reanalyses to infer their impact.In this study, regionally masked precipitation is averaged over an area within an 88 radius around each WD track point (Hunt et al. 2018b(Hunt et al. , 2019a) ) for both reanalysis and satellite precipitation products.The maximum wind speed is identified by searching for the maximum grid point values within a 58 radius (geodesic) of the tracked center (detailed discussions in sections 3b, 3c, and 3d).

2) SELECTION METHODOLOGY
The India Meteorological Department provides a challenging definition of WDs to work with.It defines them as a "cyclonic circulation/trough in the mid-to lower troposphere Unauthenticated | Downloaded 06/20/22 01:48 PM UTC or as a low-pressure area at the surface, which occur in the midlatitude westerlies and originate over the Mediterranean Sea, Caspian Sea, and the Black Sea and move eastward across north India" (http://imd.gov.in/section/nhac/termglossary.pdf).The indeterminate nature of the definition motivated us to carry out track filtering for each season, compiling the best possible western disturbance/depression catalog for the region.The selection criteria consist of the following three levels of filtering: 1) First, only those tracks that sustained themselves for at least one day and traveled a minimum distance of 1000 km are considered for the study.
2) The filtered tracks obtained at level 1 then go through the second round of filtering, which involves identifying those tracks that pass through our region of interest, defined as a lat-lon box (to be called "impact boxes" henceforth) for each of the subregions (Table 1).
3) The third and final filter is the one that takes into consideration the definition of WDs provided by the IMD.The tracks from level 2 having their genesis in and around the primary moisture source for these storms, that is, the Mediterranean Sea, the Black Sea, the Caspian Sea, as well as the Atlantic Ocean (hereafter referred to as the "genesis box"; Table 1), are separated to create the WD catalog for each subregion.Figure 2 shows the three-level filtering of tracks for the KR impact box for all the three reanalysis datasets and both the subperiods.
The level-3 filtered tracks for other subregions are depicted in Figs.S1-S3 in the online supplemental material.The study primarily uses tracks from level 2 (e.g., Fig. S4) to generate WD-associated statistics for the KH regions; unless otherwise mentioned, these include tracks that may have genesis outside the defined genesis box.This will reduce seasonal data loss, which may arise because of the tracks that substantially impact the region getting missed out due to genesis box restrictions.However, most of the significant level-2 tracks originate within the genesis box and have primary control over the WD activity of the KH, as we show in this paper.The obtained tracks were manually compared with the observed WDs of the last few decades (Dimri et al. 2015).Out of the 26 reported WDs mentioned within our study period , 22 ERA5 tracks matched the actual dates of their occurrence with some minor deviations.The remaining four are probably not found due to the filtering restrictions of origin, lifetime, and minimum distance traveled.The tracks of the matched WDs are shown in Fig. S5.

3) EADY GROWTH RATE
The maximum Eady growth rate is a measure of baroclinic instability.It follows the maximum growth rates for the configuration of the Eady problem (Eady 1949).The Eady growth rate s E (s 21 ) is given by (Lindzen and Farrell 1980;Vallis 2006;Simmonds and Lim 2009) where "fcor" is the Coriolis parameter, U(z) is the vertical profile of the eastward wind component, z is the vertical coordinate (geopotential), and N is the Brunt-Väisälä frequency given by where g is the acceleration due to gravity, and u is the potential temperature.The vertical derivatives are computed as a finite difference of values at 400-and 300-hPa levels.This allows us to investigate baroclinic instability between the two pressure levels reported to be most active for WDs (Hunt et al. 2019b).
The present study applies this expression at a 6-hourly temporal resolution for all the three reanalysis datasets.

4) DYNAMICALLY COUPLED REGIONAL CLIMATE MODEL REMO GLACIER
To assess the impact of WDs on KH glaciers, the regional mass-balance estimates for the subregions are taken from a high-resolution dynamically coupled glacier-climate model REMO glacier .The atmospheric regional model REMO (Jacob and Podzun 1997) is augmented with an online dynamical glacier parameterization scheme (DGS) (Kotlarski 2007;Kotlarski et al. 2010).DGS represents surface glacier cover on a subgrid scale and calculates the energy and mass balance of the glacierized part of a grid box.This part is allowed to grow and shrink dynamically depending on the simulated mass balance but is restricted to the respective grid box's total land surface area.The integration of REMO glacier is done over the South Asian domain (78-388N, 658-988E) for the period 1989-2016, keeping in mind the availability of glacier inventory to initialize the model and used in studies such as those of Pfeffer et al. (2014) and Frey et al. (2014).The horizontal resolution of the setup is 0.228 3 0.228 (∼25 km) with 27 vertical levels and is forced by the ERA-Interim dataset as the lateral boundary conditions.Lateral boundaries are updated every 6 h and interpolated to a 2-min time step (Kumar et al. 2015).Although the simulation has been done for the whole of South Asia, the present study uses the regionally masked mass-balance values for the subregions as defined in section 2a.Some limitations of the model are discussed in section 4.

a. Frequency and track matching of WDs
After generating the cyclone tracks for the three separate reanalysis datasets, the tracks were compared using a matching methodology.The tracks for a pair of reanalyses were said to match if at least 10% of the track points overlap in time and the mean separation between them is within a 48 radius (geodesic).Figures 3a-c show the intensity (vorticity) distribution of matching and nonmatching tracks among pairs of the three reanalysis datasets in the Northern Hemisphere.The matching results show that the tracking algorithm captures the major storms (having higher intensity values) almost identically for all three datasets.The discrepancies mainly arise for low-intensity systems whose accurate identification becomes challenging because of smaller travel distances, lifetime, and differences between the reanalysis datasets.On average, 54.78% of tracks matched among each pair of reanalyses.Further sensitivity analysis using stricter criteria revealed lower matches among the datasets.The details of the pairwise matching are provided in Table 2.
Figures 3d-g illustrate the ensemble interannual frequency of WDs over KH for the different regions.Overall, it shows that KH has had only small changes in the frequency of WDs over the full period, on par with other studies (Midhuna et al. 2020).The means for different subperiods are compared using a two-tailed unpaired t test.While the KR and CH have witnessed a slightly decreased frequency, WH and EH have observed the opposite.However, none of the ensemble trends are statistically significant at a 5% level (Tables T1-T5 in the online supplemental material).We noticed an apparent twofold difference between the frequencies of KR and WH and found two reasons for this.First, the WH impact box covers almost 60% more area than the KR impact box; thus, more tracks can satisfy the level-2 filtering.Second, the WH box lies directly in the path of traveling WDs, with a significant The large interannual variability (standard deviations in Table T1) of WD frequencies and the factors controlling their numbers are discussed in section 3f.A previous study (Hunt et al. 2018a) suggested a strong relationship between WD occurrences and the corresponding jet latitudinal position for each season.Since the frequency trend for KR is not statistically significant, the abnormal increase or stability in the mass balance of Karakoram glaciers after the year 2000 cannot be solely explained by the changes in the number of systems passing through the region.Thus, investigating the associated impact is required to explain the increase, namely the WD associated precipitation.

b. Intensity of WD associated precipitation
Despite insignificant changes in WD frequencies during P1 (Figs. 3d,e; see also Tables T2 and T3), the precipitation intensity in KR and WH is seen to show a concurrent steady decline (Figs.4c,d), which kept the precipitation volume (Figs.4a,b) and, ultimately, the mass balance in check.However, the two subregions have witnessed a revival of WDs precipitation intensity during P2.The seasonal area-averaged precipitation intensity (Figs.4c,d) derived from the precipitation values added to each point of simulated tracks (refer to section 2c) provide a clearer picture behind the sudden change in sign of mass-balance rates for the KR and nearby regions.In fact, for KR, the reversed trend in P2 is steep enough to switch the overall trend toward a positive direction, unlike WH, whose overall trend stabilized around zero.Thus, the increase in precipitation intensity presents itself as one of the possible reasons behind the KR's recent glacier surge.The precipitation trends for the Karakoram (Figs. 4a,c) also closely resemble the results shown in a previous study that compared CFSR precipitation data (dynamically downscaled using WRF) against the available Pakistani in situ   It is noticeable that apart from CH, all the other subregions show an increasing trend for snowfall in P2 (Fig. S6).The WD associated precipitation intensity of KR saw a rise of ∼10% in P2 compared to P1 (Table T2).The lower panels of Fig. 4 show the spatial difference (P2 minus P1) between the mean statistics generated using masked precipitation as the intensity parameter for KR and WH.Stippling indicates a statistically significant change at 5% level.A statistically significant increase in the precipitation intensity dominates KR, whereas more than 50% of WH shows positive values, but it is not statistically significant over the highest elevated regions most likely to have glaciers.Table 3 lists the mean precipitation intensities for the top 25th, 50th, 75th, 90th, and 95th percentiles of WDs for KR and WH.It reveals how the most intense storms of the KR region (i.e., the top 90th and 95th percentile WDs) have significantly increased their mean precipitation intensities, although the same cannot be said for WH, whose increase is marginal.In contrast, CH and EH have witnessed a decrease in mean precipitation intensities for their most intense WDs (Table T6).
The decreasing influence of WDs over CH is evident, where negative values dominate most of the region (Fig. S6, lower panels).Earlier studies (Cannon et al. 2015;Dimri et al. 2015;Krishnan et al. 2019;Norris et al. 2017) have also shown that WD-associated precipitation intensities are weakening over CH, with a subsequent decrease in precipitation.Tables T2-T5 provide the various WD-associated ensemble statistics and trends, including the mean seasonal precipitation rates across the study period for all the subregions for each year.Individual reanalysis trends and standard deviation for precipitation intensity are provided in Table T1 for each subregion, with EH precipitation intensity for MERRA-2 being the only statistically significant change.

c. Volume of WD associated precipitation
None of the earlier studies have quantified the impact of WDs on the KH glacier mass balance in recent years.However, a recent study showed that WDs contribute about 80% of the winter precipitation (DJF) in the western Himalayas (Midhuna et al. 2020).Another study has shown that winter snowfall contributes 92% and 42% of the total precipitation between December-May in the catchments of the Karakoram (Shimsal) and central Himalayas (Langtang), respectively (Bonekamp et al. 2019).The present study quantifies the volumes of solid and total WD-associated precipitation (in billions of cubic meters) for each season between 1980 and 2019.It should be noted that this section uses an ensemble of ERA5 and MERRA-2 datasets since NCEP-CFSR/CFSv2 does not provide snowfall data.For every time step a track overlaps with an impact box, the associated masked precipitation (or snowfall) within the box was iteratively added to find a system's contribution to the total precipitation (or snowfall).Figure 5 depicts the WD snowfall volumes and various contributions over KR and WH.Both subregions show an overall decreasing trend in WD snowfall volume over the last four decades, but a striking reversal in trend is apparent during P2 compared to P1 (Figs. 5a,b).The enhanced precipitation intensity due to WDs (Fig. 4) probably boosted the snowfall amounts received by the subregions.However, although an increase in solid precipitation can facilitate glaciers' accumulation, it does not necessarily result in a positive mass balance.Higher precipitation or snowfall volumes may either neutralize a high ablation during the previous summer or act as a buffer for the upcoming summer months (Kumar et al. 2019).
The WD snowfall contributes about ∼65% (Fig. 5c) to total seasonal snowfall and about ∼53% (Fig. 5e) to total seasonal precipitation in the KR impact box (Table T2).However, these numbers are expected to increase further if we consider the core regions of KR having the maximum number of glaciers.Similarly, for WH, the WD snowfall contribution is ∼65% (Fig. 5d) and ∼46% (Fig. 5f) to the total seasonal snowfall and precipitation volumes, respectively (Table T3).It is evident from these numbers that snowfall is a significant contributor to regional precipitation for KR and WH during winter, and WDs are the primary source.The CH (∼64% and ∼35%) and EH (∼60% and ∼28%) regions also receive major snowfall from WDs during the same period but significantly less than KR and WH in volume (Fig. S7 and Tables T4  and T5).
Interestingly, the WD snowfall contribution to the WD precipitation for KR and WH has remained more or less stagnant.Another interesting observation is that despite a statistically significant rise in total seasonal snowfall for KR during P2, the snowfall volume from non-WD sources has seen a simultaneous and statistically significant decline (Tables T2  and T3).The non-WD snowfall volume is found to have decreased by ∼17% in KR and by ∼15% in WH during P2 (Figs. 5g,h); also, the non-WD snowfall contributions to total seasonal snowfall reveal a decreasing trend for both KR and WH (Figs. 5i,j).The increase in WD-associated snowfall coupled with a simultaneous decrease in non-WD sources coincides with the anomaly period, suggesting a crucial role of WDs in the regional mass-balance estimates.
Another interesting aspect can be observed from Fig. 6, which shows the monthly ensemble total and WD associated precipitation/snowfall volumes for KR and WH along with the ensemble mean monthly number of WDs between 1980 and 2019.The monthly WD frequencies show a gradual and continual rise in the volumes, suggesting more intense storms in terms of precipitation/snowfall rates in the latter half of the season (Fig. 6, left panels), also indicated by Norris et al. (2015).The increase in volumes can be explained through the annual cycle of regionally masked snowfall intensity and temperature (Fig. 6, right panels).The snowfall intensity gradually increases from November to March, pointing toward an increase in snowfall magnitude without the need for increased WD activity (Norris et al. 2015).The blue bars represent the months having a higher possibility of solid precipitation since the mean temperatures are less than 08C.The annual cycle of snowfall intensity has a bimodal structure, with another peak during the later stages of the summer season, fed by ISMR (Fig. 6, right panels).The months accompanied by warmer temperatures (orange bars) are dominated by liquid precipitation, having subdued contribution toward accumulation and, subsequently, on the glaciers' mass balance.Similar patterns were observed for CH and EH as well (Fig. S8).

d. Intensity of WDs based on wind speed
The increase in WD-related impacts over KR and WH has resulted in an enhanced precipitation rate, volume, and snowfall contribution in P2.With the emergence of the anomaly over the Karakoram and no visible change in WD frequencies, it becomes important to analyze the strength of WDs (in terms of wind speed) to investigate any potential link between them.It is challenging to categorize WDs since IMD does not provide a wind speed-based category-wise distribution table for them.This section attempts to provide a benchmark for various intensity levels of the tracked WDs based on wind speed intensity.The filtered tracks based on the IMD definition of genesis (level 3) are used for classification.First, the vertically averaged wind speed of 400-300-hPa levels is added to each point of the tracks throughout its lifetime by searching for the maximum value of the full-resolution grid points within a 58 radius (geodesic) of a WD center.For each subregion, the highest sustained wind speed achieved by track points lying inside the corresponding impact boxes are used.Using simple percentiles, WDs are classified into five different intensity levels (ILs) described in Tables T7-T10.The magnitude of wind speeds for each intensity level is provided in Table T11.Figure 7 illustrates the structure of the percentage of tracks falling into different IL divisions for each of the three reanalysis datasets.All four subregions achieve their peak at IL-3, making it the most crucial level in terms of frequency and impact.A previous study has shown how the associated precipitation of extratropical cyclones can be explained using a simple scaling relationship where the predictor is the product of mean near-surface wind speed and total column water vapor (Pfahl and Sprenger 2016).It further revealed that wind speed is the dominant predictor in the regions that are typical of sufficient moisture availability.Thus, stronger winds may be associated with enhanced precipitation, directly impacting the accumulation of KH glaciers.Interestingly, the ensemble values reveal that the subregions have witnessed an increase in IL-3 and IL-4 WDs during P2 (Tables T7-T10), barring the Karakoram, whose numbers seem to be datadependent.This division hints toward a gradual increase in higher intensity WDs over the Karakoram and western Himalayas, further influencing the volume of precipitation (mainly snowfall) received by the subregions.

e. Dynamics controlling the genesis of WDs and their impact on the Karakoram-Himalayas
The increased WD precipitation intensity over KH can be associated with synoptic-scale changes such as an increase in the baroclinic instability (Tierney et al. 2018).Figure 8a shows the difference in the ensemble mean 6-hourly Eady growth rate (Hoskins and Valdes 1990) between P1 and P2 over the Eurasian region.The baroclinic instability is found to have increased over most of the genesis box, particularly over the Mediterranean Sea, the Black Sea, and the eastern Atlantic Ocean, which are the major moisture providers for cyclogenesis.Another lobe of positive values is seen residing over the KR, WH, and CH, extending across west Kunlun and Tibet.In contrast, EH, has seen a reduction in baroclinic instability.It remains to be seen whether changes in the baroclinicity arise from changes in the vertical shear (temperature gradient) or static stability and what role moisture plays in enhancing the baroclinicity (O'Gorman 2011).
The argument is given further weight by the difference plot of the track (Fig. 8b) and genesis (Fig. 8c) densities of all storm tracks of the Northern Hemisphere between the two subperiods where stippling depicts the regions of statistically significant change at 5% level.The regions of positive differences mostly coincide with the increased baroclinic instability, signaling a more conducive environment for cyclogenesis in P2.Substantial control of baroclinic instability can be gauged by the statistically significant Spearman's correlations exhibited at 5% level for the track (r s 5 0.63; p 5 0.000) and genesis (r s 5 0.48; p 5 0.002) densities throughout the last 40 years.However, the enhanced baroclinic instability and track/genesis densities seem to have a relatively lesser impact on WD frequencies over the subregions, as defined in section 2c(2) and depicted in Fig. 3.Moreover, certain stretches along the path of WDs (near the Caspian Sea) depict lesser control of baroclinicity over the track and genesis densities, revealing a diminishing impact of baroclinic instability while possibly signaling enhanced impact of other factors such as barotropic instability in modulating the densities.These possible other factors require further investigation for their complete identification.There is a strong possibility of more storms being generated in the genesis box but not reaching the subregions, thus not being reflected in the frequency analyses, although the increased baroclinic instability seems to impact WD intensities (Fig. 7) during P2.Baroclinicity can impact both the frequency and intensity of cyclones (Tierney et al. 2018) but is not the only factor impacting them.Barotropic instability can also play a part in WD development (Hunt et al. 2018a;Lee and Kim 2003), but more detailed studies are required to understand the relevant roles of both types of instabilities.Another critical factor that controls the count of extratropical cyclones is the mean latitudinal position of the subtropical westerly jet (STWJ; Hunt et al. 2018a).The eastbound STWJ is the carrier of WDs, and the mean residing latitudinal position of the STWJ determines the path taken by the WDs. Figure 9 depicts the distribution plots of the ensemble mean jet latitudes, the ensemble mean jet intensity, and the ensemble maximum jet intensity over the genesis and Himalayan boxes (Table 1) for both subperiods.Kernel density estimation (KDE) is used to compute the smoothed distributions.The winds of individual reanalyses are averaged before applying KDE.The curves have a noticeable structure, with multiple visible modalities.Even though this method smooths the curve, the multiple modalities are still retained.Figure 9a suggests that the mean jet stream latitude position over the genesis box has shifted poleward, establishing itself closer to a more moisture-abundant Mediterranean Sea, facilitating the generation of stronger WDs.
The distributions over the Himalayas, interestingly, exhibit a remarkable change.In P1, the shift of the jet is apparent with two distinct peaks (Fig. 9b).However, the jet seems to have anchored itself with more uniform probability values in P2, suggesting a broader mean jet latitudinal bandwidth within a season and an increased probability of the jet to reside at all the latitudes across the spread.Thus, unlike P1, there seems to be a strong possibility of a slower jet migration toward the north during P2.A delayed northward shift of the mean jet across northern India is shown to modulate the interannual variability of WD frequencies (Hunt et al. 2018a).The middle panels of Fig. 9 depict the distribution of mean jet intensities, showing a slight increase in jet speeds over the genesis box but a significantly higher increase over the Himalayas.In contrast, the maximum jet intensities (Fig. 9, bottom panels) show a significant decrease over the genesis box but a slight increase over the Himalayan box.However, the multiple modalities in the distributions suggest variability in propagation speeds of WDs within a season.Changes in the jet location and speed may also influence the baroclinicity via changes in the vertical shear.All the distributions presented, except Figs.9c and 9f, were found to be significantly different between P1 and P2 using the Kolmogorov-Smirnov significance test (Table T12).
To assess the increase in WD strength and moisture-holding capacity, moisture parameters such as saturated vapor pressure (svp), vapor pressure (vp), and cloud liquid water content (clwc) are analyzed over the genesis box. Figure S9 depicts the differences in moisture parameters between P1 and P2.The svp and vp are calculated using the Clausius-Clapeyron equations.All three parameters have witnessed an increase in their mean values, with positive differences dominating across the extent of the genesis box.The increase in the moisture content enhances baroclinicity and facilitates the generation of more intense and possibly a higher number of WDs (Durran and Klemp 1982).
f. Quantifying the impact of WDs on the Karakoram anomaly The KH region has a very sparse network of meteorological data recording stations, mostly placed in valleys or lower altitudes.The region's rain gauge data (or snowfall) are highly  , whereas the blue shades are for P2 .
Unauthenticated | Downloaded 06/20/22 01:48 PM UTC unreliable (Qin et al. 2009;Norris et al. 2015Norris et al. , 2017)), and therefore the observed precipitation estimates are likely underestimated, not only because of the undercatch but also because of the lack of a general factor of increase in precipitation with altitude (Kumar et al. 2015;Palazzi et al. 2013;Winiger et al. 2005;Norris et al. 2017).Dynamical regional climate models such as REMO glacier act as capable proxies for such a complex region (Kumar et al. 2015(Kumar et al. , 2019)).Coupled with a dynamical glacier scheme (DGS) (Kotlarski 2007), REMO glacier simulates the glacier mass balance on each grid box's dynamically adjusted surface fraction depending on the accumulation and ablation conditions (Kotlarski et al. 2010).
Figure 10 shows the spatial and time series plots for the REMO glacier simulated mass balance for the Karakoram region.The top panel is the annual mass balance estimates calculated for a mass balance year (October to the following September), the middle panel is for the winter season, which also corresponds with the five months considered in our study [November-March (NDJFM)], and the bottom panel is for the remaining months [April-October (AMJJASO)].The simulation period is from 1989 to 2016.It is noticeable how a positive mass balance (blue) dominates the NDJFM season (Fig. 10b), whereas negative values (red) dominate the rest of the year (the AMJJASO season; Fig. 10c), clearly behaving as expected.However, the annual spatial plot suggests that a few pockets of the southern Karakoram retain their surging/stable status (Fig. 10a).Thus, it clearly shows how the accumulation during the NDJFM season plays a vital role in the annual mass budget estimation for the region.
The present study has established that WDs are the major contributor to snowfall in the Karakoram region (∼65%), and accumulation is directly impacted by the snowfall received, which explains about ∼60% of glacier mass-balance variability (Kumar et al. 2019).The annual time series plots also suggest that the region has approached stability in the last two decades, and the NDJFM mass balance is the critical player in its enhancement, since although the trend for AMJJASO is positive in P2, the mass balance values are negative (Fig. 10, right panels).The positive trend of mass balance in AMJJASO can be attributed to a couple of factors (only the trend is positive;   and P2 (2001-19).The baroclinic instability is examined between the 300-and 400-hPa levels.Also shown are differences between P1 and P2 in (b) track density and (c) genesis density.Stippling indicates significant change at 5% level.
the values are still negative).First, the Karakoram has been less sensitive to warming in recent decades during summer (Kapnick et al. 2014;Kumar et al. 2019); second, enhanced winter accumulation increases the albedo, which is one of the most crucial factors controlling the glacial melt.Both these factors combine to impact the trend of the summer mass budget positively (Kumar et al. 2019).
As already established in the above discussion, WDs feed the accumulation, ultimately impacting the regional mass balance.Therefore, the best way to quantify the impact of WDs on glaciers is to determine the precipitation ratios by dividing the precipitation amount received due to WDs by total precipitation, masked over the glaciated region.The glacier fraction data used in the present study have been extensively used in studies such as those of Siderius et al. (2013) and Kumar et al. (2015Kumar et al. ( , 2019)).The WD-associated precipitation intensity obtained from the track statistics was used to calculate the amount of WD precipitation accumulated over the glaciated fraction.The total precipitation over the glaciated fraction was calculated next to quantify the exact contribution of WD precipitation.The resulting numbers (in percentages) reveal a statistically significant increase in the contribution ratio for the region from P1 to P2 (Fig. 10h).The ensemble mean contribution of WDs over the glaciated region rose from ∼37% in P1 to about ∼47% in P2, a relative jump of about ∼27%.The rise was found to be statistically significant at the 1% level (p value 5 0.0077).The contribution of WDs in some recent years has gone beyond 60%, clearly establishing how vital the WDs are for the recent glacial surge/stability observed in the region.The REMO glacier simulated mass balances for other subregions are shown in Fig. S10.A previous study has also conducted a spatial analysis of mass-balance estimations and showed that the model was able to capture the Karakoram anomaly reasonably well (Kumar et al. 2019).

Summary and conclusions
This study's objective is to assess the direct impact of western disturbances (WDs) on the regional winter precipitation patterns of the Karakoram-Himalayan region (KH), which control the accumulation process and, ultimately, the mass balance of its glaciers.By quantifying and analyzing WD frequency, intensity, precipitation/snowfall volumes, and wind speed, the study highlights the role played by WDs in the emergence of the "Karakoram anomaly" in recent decades.The study also analyzed the possible dynamical factors behind this enhanced WD activity over the region.To achieve the objectives, a tracking algorithm has been applied to the 3-hourly upper-tropospheric vorticity field for 39 seasons (November-March) of three separate reanalysis datasets to understand the impact of the contribution of WDs to the regional mass-balance variability of Karakoram-Himalayan glaciers for two subperiods: P1  and P2 (2001-19).Based on the definition of WDs provided by IMD, the simulated tracks underwent three-level filtering, providing a possible catalog of the synoptic-scale phenomena passing through  , whereas the solid black line represents the KDE for P2 ."Jet stream latitude" is defined as the 3-hourly latitude of maximum upper-level (averaged between 400 and 300 hPa) zonal winds zonally averaged over the genesis box and Himalayan box as defined in Table 1 (Barnes and Hartmann 2010;Woollings et al. 2010) for the November-March season from 1980 to 2019."Mean jet stream latitude" is defined as the ensemble mean of jet stream latitude calculated separately for each of the three reanalyses.Similarly, "mean jet intensity" and "maximum jet intensity" are the ensemble mean and maximum wind speed of the identified jet latitude, respectively.the four subregions of the KH, namely the Karakoram (KR), western Himalayas (WH), central Himalayas (CH), and eastern Himalayas (EH).The WD precipitation intensity for Karakoram rose by more than 10% during P2.The increasing trend coincides with the period in which the Karakoram anomaly was identified.However, the WD frequencies remained steady across the study period in all four subregions.Despite a static trend in frequencies, the regional mass balance of KR and WH in P1 observed a downward trend due to a simultaneous decrease in WD-associated precipitation intensities.However, the trends reversed in P2 for both the subregions, with spatial analysis revealing that KR and WH are dominated by positive differences (P2 minus P1) in mean WD-associated precipitation intensities.
The present study provides a holistic view of the WD associated activities for the four major regions of the Karakoram-Himalayas, which has not been attempted before.For the first time, the contribution of WDs to the regional precipitation regime is quantified.Snowfall from WDs constituted a major share of total seasonal precipitation over the KR (∼53%) and WH (∼46%) regions, and to a lesser extent over the CH (∼35%) and EH (∼28%).The WD-associated total precipitation and snowfall volumes (in billions of cubic meters) were also computed for each subregion to gauge their impact on regional mass-balance estimations.It was established that despite an insignificant change in WD frequencies, the enhanced precipitation/snowfall intensity boosted the snowfall amount received by the Karakoram and western Himalayas in P2.The volume of WD associated snowfall in KR and WH declined steadily during P1, but a striking reversal in trends is apparent for the Karakoram and western Himalayas since the start of P2.It was found that the percentage contribution of WD snowfall to total seasonal snowfall and total seasonal precipitation also increases for KR and WH, with a simultaneous statistically significant decrease in contributions from non-WD sources.Moreover, despite the low variance in mean monthly WD frequencies, the snowfall intensity and the associated precipitation volumes are found to intensify gradually across the five months (November-March) for all the subregions.
For the first time, the simulated WDs are identified and distributed in percentile-based intensity levels (ILs) computed using their wind speed strength.The WDs were distributed into five different intensity levels of increasing wind speeds.IL-3 is the most prominent in terms of frequency and precipitation amount received by the subregions.The Karakoram and western Himalayas observed an increase of systems in IL-3 and IL-4, thus directly impacting the accumulation period's snowfall amount.The recently enhanced WD intensity over the Karakoram-Himalayas was traced back to the genesis region of WDs.On carefully analyzing cyclogenesis factors over the Mediterranean Sea using various moisture-related parameters, it was found that an increase in the moisture content of the genesis box may have impacted the intensity and duration of WDs, if not their frequencies.At the same time, the role of baroclinic instability and the mean latitudinal position of the subtropical westerly jet (STWJ) was also found to be important for WD genesis and strength.Significant changes in both the factors were observed in P2 for both the genesis and impact boxes.
A mass-balance calculation for the Karakoram-Himalayan glaciers using REMO glacier captured the Karakoram anomaly.The present study claims that the revival of WD intensity and precipitation volumes in P2 to be one of the most important drivers behind its establishment.The anomaly can be linked to the variability in precipitation patterns observed during its major accumulation period (i.e., the winters).The KR and WH witnessed a concurrent rise in WD precipitation intensity, precipitation volume, snowfall volume, percentage of snowfall contributed by WDs, and moisture carrying capacity.A recent study highlighted snowfall variability as the primary controlling factor of mass balance in the region, explaining ∼60% of the variability (Kumar et al. 2019).However, temperature variability also plays an essential role (Forsythe et al. 2017;Kapnick et al. 2014), which is not discussed in the study.Last, to quantify the direct impact of WDs on the glaciers of the Karakoram, the annual precipitation ratios were calculated by dividing the amount of WD-associated precipitation by total precipitation, masked over the glaciated fraction.The contributions jumped from ∼37% in P1 to about ∼47% in P2, a statistically significant rise of about ∼27% in impact.A strong accumulation season can impact regional mass balance in two ways.It can either nullify the effect of strong ablation during the previous summer or create a buffer for strong ablation in the upcoming summer months (Kumar et al. 2019).In both cases, a higher accumulation raises the probability of a positively skewed mass balance.In contrast, the central and eastern Himalayas continue to lose mass, although at slightly lower rates than before.
However, several questions remain unanswered, such as the trigger behind changes in dynamical and thermodynamic parameters controlling the WD genesis.Some studies have mentioned the possible eastward shift of the Karakoram anomaly toward the western Kunlun and Pamir regions (Brun et al. 2017;Berthier and Brun 2019;Gardelle et al. 2013).It remains to be seen whether the changes in the WD core genesis zone has a potential link with the shift.Another aspect would be to quantify uncertainties in climate trends derived from various reanalysis products, arising due to changes in the observing systems (Bengtsson 2004).While there may be uncertainties in snowfall from reanalyses (Orsolini et al. 2019), they provide homogeneous data with some confidence for analyzing regions such as the KH region, where there is an acute shortage of reliable and long-term observational data records (Kumar et al. 2015).For example, the precipitation data of ERA5, when compared with other precipitation products over the study region such as TRMM, CRU, and GPM, were found to systematically overestimate the seasonal precipitation (November-March) for the common period 2002-18 (Fig. S11, left panels) for all the subregions.However, these observational products also have uncertainties over high orography (Bharti and Singh 2015).Interestingly, the percentage contribution of each season to the total annual precipitation volume is in general agreement with the other datasets except for the eastern Himalayas (the bias is systematic here as well), showing that the reanalysis can capture the climatic signals on par with other observational and satellite data products (Fig. S11, right panels).
Another way to look into the performance of reanalysis datasets with respect to the available observational satellite products (e.g., TRMM) is to associate the satellite precipitation with the respective tracks of different reanalyses in the same as for the reanalysis precipitation and compare their respective intensities.,j) underestimates precipitation when compared to TRMM.However, ERA5 demonstrates a very good spatial and statistically significant correlation with TRMM (Fig. S12c), followed by MERRA-2 (Fig. S12g) and the least by NCEP-CFSR (Fig. S12k).This clearly shows the uncertainty between reanalysis datasets, in particular in regions of complex topography where in situ measurements are extremely scarce.They also have difficulty capturing the interannual variability (Figs. S12d,h,l), with ERA5 being the best among all with a systematic bias.The performance of ERA5 does suggest improvements in results when a different model formulation and assimilation technique is employed.Despite obvious shortcomings, reanalysis datasets can be considered comparable to the satellite products, bearing in mind the uncertainties in both, since our study period coincides with the satellite era, which contributes the largest number of observations.
Similarly, the simulated mass-balance estimation of REMO glacier is also subject to some apparent limitations such as an overestimation of orographic precipitation in some grid boxes due to sharp orographic gradients among adjacent grid boxes, neglect of subgrid variability of atmospheric forcing parameters for the embedded glacier scheme, and uncertainties of the observational glacier inventory used to initialize glaciers in the RCM (Kotlarski et al. 2010;Kumar et al. 2015).However, the model captures the spatial patterns of positive mass balance in some parts of the region, especially over the Karakoram (Kumar et al. 2015(Kumar et al. , 2019)).The simulated results agree with previous studies that provide confidence in the model setup and its general applicability for future climate and glacier change studies over the KH region.The uncertainty around the tracking scheme can also be explored to check for the sensitivity of the results.Using 3-hourly datasets can improve some aspects of the tracking, such as false connections, but it introduces other issues such as more frequent multiple centers.However, this is reduced to some extent by the spectral filtering, as done in the present study.Any vorticity-based scheme is likely to produce similar results as long as the identification criteria are the same.Alternatively, other fields such as upper-level geopotential or streamfunction could also be used for tracking in future work to assess the sensitivity to the identification variable.Despite some limitations and uncertainties, it should be kept in mind that the tracking results are consistent with the changes in other largescale factors, thus providing confidence in the results over the study region.The study of future WDs using scenario-based model simulations is another aspect to work on, keeping in mind the importance of glaciers and the impending future water shortage.

FIG 1 .
FIG 1.(a) The four regions defined in Bolch et al. (2012) and used in the present study.Yellow circles denote the location of available IMD stations, whereas green triangles denote the CRU land stations.The inset map shows the location and boundaries of the KH.(b) The REMO glacier simulated mass balance (in meters of water equivalent) of the KH region from 1989 to 2016, adapted from Kumar et al. (2019).The red oval highlights the anomalous positive mass balance captured by the model.
FIG 2. The WD tracks for KR impact box (in cyan).The magenta box denotes the genesis box as defined by IMD.All the tracks passing through KR having their origin in the genesis box during P1 (1980-2000; red dots) and P2 (2001-19; yellow dots) are shown for (a) ERA5, (b) MERRA-2, and (c) NCEP-CFSR/CFSv2.
FIG 3. (a)-(c) Paired intensity distribution of matching and nonmatching tracks among the three reanalysis datasets over the Northern Hemisphere.(d)-(g) Ensemble frequency time series (lines) of tracks passing through KR, WH, CH, and EH, respectively, with the 95% confidence interval (shading).Ensemble frequency is defined by averaging the frequency obtained from all three reanalyses.This helps remove multiple occurrences of the tracks that matched while still giving some weightage to those that qualify the selection criteria but did not match.The dotted lines depict the period-wise trends.
observational data records[Figs.5a and 6a of Norris et al.  (2018)].Their overall trend (1979Their overall trend ( -2010) )  for winter precipitation is not significant, similar to what our calculations also suggest.Moreover, visual analysis suggests that their trend for overlapping years of P1 (P2) is decreasing (increasing), further validating our findings.Our results also seem to reasonably capture the spatially depicted precipitation values shown in studies such as those ofCannon et al. (2015) andKrishnan et al. (2018) over regions closer to our study domain.

FIG 4 .
FIG 4. (a),(b) Ensemble-mean precipitation volume [in billions of cubic meters (bcm)] shown for the KR and WH, respectively.(c),(d) Ensemble-mean precipitation intensity (in mm day 21 ) for the KR and WH.The black line denotes the overall trend.The red dashed line denotes the trend for P1 (1980-2000), whereas the blue dashed line denotes the trend for P2 (2001-19).(e),(f) Spatial plots for the difference in mean precipitation intensity between P1 and P2, for the KR and WH.The red lines denote the extent of the impact box.Stippling indicates significant change at 5% level.

FIG 5 .
FIG 5. (a),(b) Volume of WD snowfall (in bcm) for the KR and WH, respectively.(c),(d) WD snowfall contribution (in %) in total seasonal snowfall in the KR and WH.(e),(f) WD snowfall contribution (in %) in total seasonal precipitation for the KR and WH.(g),(h) Non-WD snowfall volume (in bcm) in the KR and WH.(i),(j) Non-WD snowfall contribution (in %) in total seasonal snowfall in the KR and WH.The black line denotes the overall trend.The red dashed line denotes the trend for P1 (1980-2000), whereas the blue dashed line denotes the trend for P2 (2001-19).

FIG 6 .
FIG 6. (a),(c) Ensemble-mean monthly (for November-March) total and WD associated total precipitation volume, snowfall volume, and the mean number of tracks for 39 seasons for the KR and WH, respectively.(b),(d) Annual cycle of ensemble snowfall intensity (in mm day 21 ) and temperature for the KR and WH.

FIG 7 .
FIG 7. Percentage of tracks falling in different intensity levels for all three datasets for the (a) KR, (b) WH, (c) CH, and (d) EH regions.The brown shades are for P1, whereas the blue shades are for P2.
FIG 9. (a),(b) Mean jet stream latitudes with KDE for the genesis and Himalayan boxes.(c),(d) Mean jet stream intensity with KDE for the genesis and Himalayan boxes.(e),(f) Maximum jet stream intensity with KDE for genesis and Himalayan boxes.The orange and blue step curves depict the histograms of the distribution (bins 5 20).The dashed black line represents the KDE for P1, whereas the solid black line represents the KDE for P2."Jet stream latitude" is defined as the 3-hourly latitude of maximum upper-level (averaged between 400 and 300 hPa) zonal winds zonally averaged over the genesis box and Himalayan box as defined in Table 1(Barnes and Hartmann 2010; Woollings et al. 2010)  for the November-March season from 1980 to 2019."Mean jet stream latitude" is defined as the ensemble mean of jet stream latitude calculated separately for each of the three reanalyses.Similarly, "mean jet intensity" and "maximum jet intensity" are the ensemble mean and maximum wind speed of the identified jet latitude, respectively.

FIG 10 .
FIG 10.Spatial plot of REMO glacier simulated mass balance of Karakoram for the (a) annual (October-September), (b) NDJFM (November-March), (c) AMJJASO (April-October) periods for 1989-2016.Also shown are time series of the REMO glacier simulated mass balance for the (d) annual, (e) NDJFM, and (f) AMJJASO periods.Orange dashed lines with crosses depict the trend for P1 (1990-2000), whereas green dashed lines with diamonds depict the trend for P2 (2001-16).(h) The percentage contribution of WDs in total accumulation over the glaciated fraction of the Karakoram using the ensemble of three reanalysis datasets for the period 1980-2019.The red dashed line depicts the trend for P1 (1980-2000), whereas the blue dashed line depicts the trend for P2 (2001-19).Note that (d) is adapted from Fig. S4(c) of Kumar et al. (2019).
Figure S12 (top row) shows the spatial Unauthenticated | Downloaded 06/20/22 01:48 PM UTC plots of Karakoram WD precipitation intensity for ERA5 and TRMM (added to ERA5 tracks) along with their spatial correlation and time series.The middle panels are for MERRA-2, while the lower panels are for NCEP-CFSR.The figures suggest an overestimation of precipitation by ERA5 (Figs.S12a,b) and MERRA-2 (Figs.S12e,f), whereas NCEP-CFSR (Figs.S12i

TABLE 1 .
The table shows the lat-lon extent of all the impact boxes, the genesis box, and the Himalayan box.

TABLE 2 .
Sensitivity analysis of various matching criteria used for the generated tracks of ERA5, MERRA-2, and NCEP reanalysis datasets for Northern Hemisphere tracks from 1980 to 2019.All numbers are in %.Ensemble values are highlighted in boldface.

TABLE 3 .
Mean precipitation intensity of WDs for respective percentiles during P1 and P2 along with their difference (P2 minus P1) for KR and WH.Units are mm day 21 .Differences are set in italic for emphasis.